rm(list = ls())
cat("\014")
library(tidyverse, quietly = T)
## Warning: package 'tidyverse' was built under R version 3.5.1
## -- Attaching packages ---------------------------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0 v purrr 0.2.5
## v tibble 2.0.1 v dplyr 0.7.8
## v tidyr 0.8.1 v stringr 1.3.1
## v readr 1.3.1 v forcats 0.3.0
## Warning: package 'ggplot2' was built under R version 3.5.2
## Warning: package 'tibble' was built under R version 3.5.2
## Warning: package 'tidyr' was built under R version 3.5.1
## Warning: package 'readr' was built under R version 3.5.2
## Warning: package 'purrr' was built under R version 3.5.1
## Warning: package 'dplyr' was built under R version 3.5.2
## Warning: package 'stringr' was built under R version 3.5.1
## Warning: package 'forcats' was built under R version 3.5.1
## -- Conflicts ------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(sf, quietly = T)
## Warning: package 'sf' was built under R version 3.5.2
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(gganimate, quietly = T)
library(openair, quietly = T)
## Warning: package 'openair' was built under R version 3.5.1
library(lubridate, quietly = T)
## Warning: package 'lubridate' was built under R version 3.5.1
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
library(zoo, quietly = T)
## Warning: package 'zoo' was built under R version 3.5.1
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(ggthemes, quietly = T)
## Warning: package 'ggthemes' was built under R version 3.5.2
library(ggmap, quietly = T)
## Warning: package 'ggmap' was built under R version 3.5.2
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(sp, quietly = T)
## Warning: package 'sp' was built under R version 3.5.1
library(magick, quietly = T)
## Warning: package 'magick' was built under R version 3.5.2
## Linking to ImageMagick 6.9.9.14
## Enabled features: cairo, freetype, fftw, ghostscript, lcms, pango, rsvg, webp
## Disabled features: fontconfig, x11
Import the time-activity diary for the line and tidy it up.
diary <- read_csv('air_quality_data/victoria_line_2014-12-08.csv', col_types = cols())
base <- tibble(date_time = seq(min(diary$start_time), max(diary$end_time), by = 'min'))
base <- left_join(base, diary, by = c('date_time' = 'start_time')) %>%
select(-fake_id)
base <- left_join(select(base, -end_time),
filter(base, date_time != end_time),
by = c('date_time' = 'end_time')) %>%
mutate(line.x = if_else(is.na(line.x) & !is.na(line.y), line.y, line.x),
station.x = if_else(is.na(station.x) & !is.na(station.y), station.y, station.x),
environment.x = if_else(is.na(environment.x) & !is.na(environment.y), environment.y, environment.x),
sub_environment.x = if_else(is.na(sub_environment.x) & !is.na(sub_environment.y), sub_environment.y, sub_environment.x)) %>%
select(date_time, line.x, station.x, environment.x, sub_environment.x) %>%
rename(line = line.x, station = station.x, environment = environment.x, sub_environment = sub_environment.x)
## Warning: package 'bindrcpp' was built under R version 3.5.1
rm(diary)
Import the measured air quality to join to the time-activity diary
air_quality <- read_csv('air_quality_data/LUexportMay2016.csv', col_types = cols()) %>%
rename_all(tolower) %>%
mutate(datetime = as.POSIXct(datetime, format='%d/%m/%Y %H:%M')) %>%
filter(sitecode == 'CAR')
Join the time activity and air quality
line <- left_join(base, air_quality, by = c('date_time' = 'datetime')) %>%
mutate(scaledvalue = if_else(species == 'PM25', value, scaledvalue)) %>%
select(-value) %>%
rename(concentration = scaledvalue)
rm(air_quality, base)
Import the background air quality for that day to correct the PM2.5 measurements
background_pm25 <- importKCL(site = "kc1", year = c(2014,2015,2016), pollutant = "pm25", met = FALSE,
units = "mass", extra = FALSE)
## NOTE - mass units are used
## ug/m3 for NOx, NO2, SO2, O3; mg/m3 for CO
## PM10_raw is raw data multiplied by 1.3
background_pm25 <- data.frame(background_pm25, day = as.Date(format(background_pm25$date)))
background_pm25 <- aggregate(pm25 ~ day, background_pm25, mean)
Join the background PM2.5 to the monitored concentrations
line <- mutate(line, day = date(date_time)) %>%
left_join(background_pm25, by = c('day' = 'day')) %>%
select(-day)
Correct the PM2.5 data
line <- mutate(line, concentration =
case_when(concentration > pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ pm25 + (1.82 * (concentration - (pm25/0.44))),
concentration <= pm25/0.44 & species == 'PM25' & !is.na(concentration) ~ concentration * 0.44,
TRUE ~ concentration)) %>%
select(-pm25)
Get the station locations
stations <- st_read('https://raw.githubusercontent.com/dracos/underground-live-map/master/bin/stations.kml', quiet = T) %>%
select(-Description) %>%
rename_all(tolower) %>%
mutate(name = gsub(' Station', '', name)) %>%
mutate(name = gsub("'", '', name)) %>%
st_zm(drop = TRUE)
Join locations to the concentrations and station names
line <-left_join(line, stations, by = c('station' = 'name')) %>% st_as_sf()
Fill in the missing coordinates using linear interpolation
st_geometry(line) <- as_tibble(st_coordinates(line)) %>%
mutate(X = na.approx(X), Y = na.approx(Y)) %>%
st_as_sf(coords = c('X', 'Y')) %>%
st_geometry()
Start trying to animate. First the graph animation.
text_placement <- filter(line,species == 'PM25') %>%
as_tibble() %>%
select(-geometry) %>%
group_by(species) %>%
summarise(max_x = max(date_time), max_y = max(concentration))
line <- left_join(line, text_placement, by = c('species' = 'species'))
graph <- ggplot(data = filter(line,species == 'PM25')) +
geom_path(aes(date_time, concentration), colour = '#0099CC', size = 1.2) +
geom_point(aes(date_time, concentration), colour = 'red', size = 4) +
geom_text(aes(max_x, max_y,label = round(concentration,0)),
size = 10, hjust = 1, vjust = 1, colour = '#0099CC') +
theme(axis.text.x = element_blank(),
axis.text.y = element_text(size = 12),
axis.title.y = element_text(size = 12),
axis.line = element_line(colour = 'black')) +
xlab('') +
ylab('PM2.5 ug/m3') +
transition_reveal(date_time)
graph_animation <- animate(graph, fps = 3, height = 600)
anim_save(file='graph_animation.gif', animation = graph_animation, path = 'outputs')
Now the spatial animation.
line <- st_set_crs(line,4326)
map_area <- c(as.vector(st_bbox(line))[1]-0.01, as.vector(st_bbox(line))[2]-0.01, as.vector(st_bbox(line))[3]+0.01, as.vector(st_bbox(line))[4]+0.01)
map <- ggmap(get_map(location = map_area), source = "google", maptype = "roadmap", zoom = 13) +
geom_path(data = as_tibble(st_coordinates(filter(line,species == 'PM25')[1:56,])), aes(X, Y), colour = "#0099CC", size = 1) +
geom_sf(data = filter(line,species == 'PM25'), size=4, colour = 'red',inherit.aes = FALSE) +
coord_sf(datum=NA) +
transition_time(date_time) +
ease_aes('linear') +
theme(axis.title = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(colour = 'black', fill = NA),
axis.text = element_blank())
## Source : http://tile.stamen.com/terrain/13/4092/2720.png
## Source : http://tile.stamen.com/terrain/13/4093/2720.png
## Source : http://tile.stamen.com/terrain/13/4094/2720.png
## Source : http://tile.stamen.com/terrain/13/4092/2721.png
## Source : http://tile.stamen.com/terrain/13/4093/2721.png
## Source : http://tile.stamen.com/terrain/13/4094/2721.png
## Source : http://tile.stamen.com/terrain/13/4092/2722.png
## Source : http://tile.stamen.com/terrain/13/4093/2722.png
## Source : http://tile.stamen.com/terrain/13/4094/2722.png
## Source : http://tile.stamen.com/terrain/13/4092/2723.png
## Source : http://tile.stamen.com/terrain/13/4093/2723.png
## Source : http://tile.stamen.com/terrain/13/4094/2723.png
## Source : http://tile.stamen.com/terrain/13/4092/2724.png
## Source : http://tile.stamen.com/terrain/13/4093/2724.png
## Source : http://tile.stamen.com/terrain/13/4094/2724.png
## Source : http://tile.stamen.com/terrain/13/4092/2725.png
## Source : http://tile.stamen.com/terrain/13/4093/2725.png
## Source : http://tile.stamen.com/terrain/13/4094/2725.png
## Source : http://tile.stamen.com/terrain/13/4092/2726.png
## Source : http://tile.stamen.com/terrain/13/4093/2726.png
## Source : http://tile.stamen.com/terrain/13/4094/2726.png
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
map_animation <- animate(map, fps = 3, height = 600)
anim_save(file='map_animation.gif', animation = map_animation, path = 'outputs')
graph_animation
map_animation